// File:    func.cpp
// Author:  Brian Allen Vanderburg II
// Purpose: ExprEval internal functions
//------------------------------------------------------------------------------

// Includes
#include <new>
#include <memory>

#include <cmath>
#include <ctime>

#include "COMMON/expr/defs.h"
#include "COMMON/expr/funclist.h"
#include "COMMON/expr/node.h"
#include "COMMON/expr/error.h"

namespace ExprEval
{

// Absolute value
//--------------------------------------------------------------------------
class abs_FunctionNode : public FunctionNode
{
public:
    abs_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        return fabs (m_nodes[0]->Evaluate (errors));
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        if(root.first>=0.)
        {
            return root;
        }
        else
        {
            return std::make_pair (-root.first, -root.second);
        }
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("not tested");
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        if(root.first>=0.)
        {
            return root;
        }
        else
        {
            return std::make_pair (-root.first, -root.second);
        }
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("abs (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};


class abs_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "abs";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new abs_FunctionNode (expr);
    }
};

// Modulus
//--------------------------------------------------------------------------
class mod_FunctionNode : public FunctionNode
{
public:
    mod_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double result;

        // Check for math errors
        errno = 0;

        result = std::fmod (m_nodes[0]->Evaluate (errors), m_nodes[1]->Evaluate (errors));
        if(errno)
        {
            errors.add_error(errno, "mod");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'mod' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'mod' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("mod (" + m_nodes[0]->get_renamed_expression (var_map_name) + ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class mod_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "mod";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new mod_FunctionNode (expr);
    }
};

// Integer part
//--------------------------------------------------------------------------
class ipart_FunctionNode : public FunctionNode
{
public:
    ipart_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double result;

        modf(m_nodes[0]->Evaluate (errors), &result);

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'ipart' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'ipart' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("ipart (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class ipart_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "ipart";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new ipart_FunctionNode (expr);
    }
};

// Fraction part
//--------------------------------------------------------------------------
class fpart_FunctionNode : public FunctionNode
{
public:
    fpart_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double dummy;

        return modf(m_nodes[0]->Evaluate (errors), &dummy);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'fpart' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'fpart' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("fpart (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class fpart_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "fpart";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new fpart_FunctionNode (expr);
    }
};

// Minimum
//--------------------------------------------------------------------------
class min_FunctionNode : public FunctionNode
{
public:
    min_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, -1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        std::vector<Node*>::size_type pos;

        double result = m_nodes[0]->Evaluate (errors);

        for (pos = 1; pos < m_nodes.size(); pos++)
        {
            double tmp = m_nodes[pos]->Evaluate (errors);
            if(tmp < result)
                result = tmp;
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'min' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'min' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        std::ostringstream os;
        os << "min (" << m_nodes[0]->get_renamed_expression (var_map_name) << ", ";

        for (std::size_t pos = 1; pos < m_nodes.size(); pos++)
        {
            os << m_nodes[pos]->get_renamed_expression (var_map_name);
            if (pos != m_nodes.size() - 1)
            {
                os << ", ";
            }
        }

        os << ")";
        return os.str();
    }
};

class min_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "min";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new min_FunctionNode (expr);
    }
};

// Maximum
//--------------------------------------------------------------------------
class max_FunctionNode : public FunctionNode
{
public:
    max_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, -1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        std::vector<Node*>::size_type pos;

        double result = m_nodes[0]->Evaluate (errors);

        for (pos = 1; pos < m_nodes.size(); pos++)
        {
            double tmp = m_nodes[pos]->Evaluate (errors);
            if(tmp > result)
                result = tmp;
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'max' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'max' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        std::ostringstream os;
        os << "max (" << m_nodes[0]->get_renamed_expression (var_map_name) << ", ";

        for (std::size_t pos = 1; pos < m_nodes.size(); pos++)
        {
            os << m_nodes[pos]->get_renamed_expression (var_map_name);
            if (pos != m_nodes.size() - 1)
            {
                os << ", ";
            }
        }

        os << ")";
        return os.str();
    }
};

class max_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "max";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new max_FunctionNode (expr);
    }
};

// Square root
//--------------------------------------------------------------------------
class sqrt_FunctionNode : public FunctionNode
{
public:
    sqrt_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        const double result = sqrt (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "sqrt");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;
        const std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);

        const double value = sqrt (root.first);
        if(errno)
        {
            errors.add_error(errno, "sqrt");
        }

        /*
         * h=sqrt (f) ==> dh/dx = (1/2*sqrt (f))df/dx
         */
        double dvalue = root.second/(2.*value+TOLERANCE);
        if(errno)
        {
            errors.add_error(errno, "sqrt");
        }
        if(errno == ERANGE)
        {
            dvalue = INF_VALUE;
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=sqrt (f) ==> d^2h/dxdy = d ((1/2sqrt (f))df/dy)/dx
         * = -(1/4)*f^(-3/2)df/dx*df/dy+(1/2)*f^(-1/2)d^2f/dxdy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double f_1  = 1/sqrt(f);
        const double f_3  = f_1 / f;
        const double df   = -0.25 * f_3 * f_x * f_y + 0.5 * f_1 * f_xy;
        return std::make_pair (sqrt (f), df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("sqrt (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class sqrt_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "sqrt";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new sqrt_FunctionNode (expr);
    }
};

// Sin
//--------------------------------------------------------------------------
class sin_FunctionNode : public FunctionNode
{
public:
    sin_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = sin (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "sin");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        if(errno)
        {
            errors.add_error(errno, "sin");
        }

        double value = sin (root.first);
        if(errno)
        {
            errors.add_error(errno, "sin");
        }

        /*
         * h=sin (f) ==> dh/dx=cos (f)df/dx
         */
        double dvalue = cos (root.first)*root.second;
        if(errno)
        {
            errors.add_error(errno, "sin");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=sin (f) ==>d^2h/dxdy = d (cos (f)df/dy)/dx
         * = cos (f)d^2f/dxdy-sin (f)df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double sinf = std::sin (f);
        if(errno)
        {
            errors.add_error(errno, "sin");
        }

        const double cosf = std::cos (f);
        if(errno)
        {
            errors.add_error(errno, "sin");
        }

        const double df   = cosf*f_xy-sinf*f_x*f_y;

        return std::make_pair (sinf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("sin (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class sin_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "sin";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new sin_FunctionNode (expr);
    }
};

// Cosine
//--------------------------------------------------------------------------
class cos_FunctionNode : public FunctionNode
{
public:
    cos_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = cos (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        double value = cos (root.first);
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        double dvalue = -sin (root.first)*root.second;
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=cos (f) ==>d^2h/dxdy = d (-sin (f)df/dy)/dx
         * = -sin (f)d^2f/dxdy-cos (f)df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double sinf = std::sin (f);
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        const double cosf = std::cos (f);
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        const double df   = -sinf*f_xy-cosf*f_x*f_y;

        return std::make_pair (cosf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("cos (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class cos_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "cos";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new cos_FunctionNode (expr);
    }
};

// Tangent
//--------------------------------------------------------------------------
class tan_FunctionNode : public FunctionNode
{
public:
    tan_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = std::tan (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "tan");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        if(errno)
        {
            errors.add_error(errno, "tan");
        }

        double cos_value = std::cos (root.first);
        if(errno)
        {
            errors.add_error(errno, "tan");
        }

        double value = std::tan (root.first);
        if(errno)
        {
            errors.add_error(errno, "tan");
        }

        double dvalue = root.second/(cos_value*cos_value+TOLERANCE);
        if(errno)
        {
            errors.add_error(errno, "tan");
        }

        if(errno==ERANGE)
        {
            if(dvalue>0)
            {
                dvalue = INF_VALUE;
            }
            else
            {
                dvalue = -INF_VALUE;
            }
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
          * h=tan (f) ==>d^2h/dxdy = d ((1/cos^2(f))df/dy)/dx
          * = 1/cos^2(f)*d^2f/dxdy+2sin (f)/cos^3(f)df/dx*df/dy
          */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double sinf = std::sin (f);
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        const double cosf = std::cos (f);
        if(errno)
        {
            errors.add_error(errno, "cos");
        }

        const double cos1f = 1./(cosf+sign (cosf)*TOLERANCE);
        const double cos2f = cos1f*cos1f;
        const double df = cos2f*f_xy+2.*sinf*cos2f*cos1f*f_x*f_y;
        return std::make_pair (sinf*cos1f, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("tan (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class tan_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "tan";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new tan_FunctionNode (expr);
    }
};

// Hyperbolic Sine
// sinh(x)=(exp(x)-exp(-x))/2
//--------------------------------------------------------------------------
class sinh_FunctionNode : public FunctionNode
{
public:
    sinh_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = std::sinh(m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "sinh");
        }
        if(errno==ERANGE)
        {
            result = INF_VALUE*sign (result);
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=sinh(f) ==> df/dx = cosh(f)df/dx
         */
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = std::sinh(root.first);
        if(errno)
        {
            errors.add_error(errno, "sinh");
        }

        double dvalue = root.second*std::cosh(root.first);

        if(errno)
        {
            errors.add_error(errno, "sinh");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
          * h=sinh(f) ==>d^2h/dxdy = d (cosh(f)df/dy)/dx
          * = cosh(f)*d^2f/dxdy + sinh(f)df/dx*df/dy
          */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double sinhf = std::sinh(f);
        if(errno)
        {
            errors.add_error(errno, "sinh");
        }

        const double coshf = std::cosh(f);
        if(errno)
        {
            errors.add_error(errno, "sinh");
        }

        const double df = coshf*f_xy+sinhf*f_x*f_y;
        return std::make_pair (sinhf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("sinh (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class sinh_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "sinh";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new sinh_FunctionNode (expr);
    }
};

// Hyperbolic Cosine
// cosh(x)=(exp(x)+exp(-x))/2
//--------------------------------------------------------------------------
class cosh_FunctionNode : public FunctionNode
{
public:
    cosh_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = std::cosh(m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "cosh");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = std::cosh(root.first);
        double dvalue = root.second*std::sinh(root.first);

        if(errno)
        {
            errors.add_error(errno, "cosh");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
          * h=cosh(f) ==>d^2h/dxdy = d (sinh(f)df/dy)/dx
          * = sinh(f)*d^2f/dxdy + cosh(f)df/dx*df/dy
          */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double sinhf = std::sinh(f);
        if(errno)
        {
            errors.add_error(errno, "cosh");
        }

        const double coshf = std::cosh(f);
        if(errno)
        {
            errors.add_error(errno, "cosh");
        }

        const double df = sinhf*f_xy+coshf*f_x*f_y;
        return std::make_pair (coshf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("cosh (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class cosh_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "cosh";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new cosh_FunctionNode (expr);
    }
};

// Hyperbolic Tangent
//tanh(x)=sinh(x)/cosh(x)
//--------------------------------------------------------------------------
class tanh_FunctionNode : public FunctionNode
{
public:
    tanh_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = std::tanh(m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "tanh");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = tanh(root.first);
        double cosh_value = cosh(root.first);
        double dvalue = root.second/(cosh_value*cosh_value);

        if(errno)
        {
            errors.add_error(errno, "tanh");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=tanh(f) ==>d^2h/dxdy = d (1/cosh^2(f)df/dy)/dx
         * = 1/cosh^2(f)*d^2f/dxdy - 2*sinh(f)/cosh^3(f)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double sinhf = std::sinh(f);
        if(errno)
        {
            errors.add_error(errno, "tanh");
        }
        const double coshf = std::cosh(f);
        if(errno)
        {
            errors.add_error(errno, "tanh");
        }

        const double cosh_1f = 1./(coshf+sign (coshf)*TOLERANCE);
        const double cosh_2f = cosh_1f*cosh_1f;
        const double df      = cosh_2f*f_xy-2.*sinhf*cosh_1f*cosh_2f*f_x*f_y;
        return std::make_pair (sinhf/coshf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("tanh (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class tanh_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "tanh";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new tanh_FunctionNode (expr);
    }
};

// Arc Sine
//--------------------------------------------------------------------------
class asin_FunctionNode : public FunctionNode
{
public:
    asin_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = asin (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "asin");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h = asin (f) ==> dh/dx = 1/sqrt (1-f^2)*df/dx
         */
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = std::asin (root.first);
        double dvalue = root.second/(std::sqrt (std::abs (1.-root.first*root.first))+TOLERANCE);

        if(errno)
        {
            errors.add_error(errno, "asin");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=asin (f) ==>d^2h/dxdy = d (1/sqrt (1-f^2)df/dy)/dx
         * = 1/sqrt (1-f^2)*d^2f/dxdy + f/(1-f^2)^(3/2)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f     = root_2.first;
        const double f_x   = root_11.second;
        const double f_y   = root_12.second;
        const double f_xy  = root_2.second;
        const double sqrtf = 1./std::sqrt (std::abs (1.-f*f)+TOLERANCE);
        if(errno)
        {
            errors.add_error(errno, "asin");
        }

        const double sqrt3f = sqrtf/(std::abs (1.-f*f)+TOLERANCE);
        if(errno)
        {
            errors.add_error(errno, "asin");
        }

        const double df      = sqrtf*f_xy+f*sqrt3f*f_x*f_y;
        return std::make_pair (std::asin (f), df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("asin (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class asin_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "asin";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new asin_FunctionNode (expr);
    }
};

// Arc Cosine
//--------------------------------------------------------------------------
class acos_FunctionNode : public FunctionNode
{
public:
    acos_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = acos (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "acos");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=acos (f) ==> dh/dx = -1/sqrt (1-f^2)df/dx
         */
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = acos (root.first);
        double dvalue = -root.second/(std::sqrt (std::abs (1.-root.first*root.first)+TOLERANCE));

        if(errno)
        {
            errors.add_error(errno, "acos");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=acos (f) ==>d^2h/dxdy = d (-1/sqrt (1-f^2)df/dy)/dx
         * = -1/sqrt (1-f^2)*d^2f/dxdy - f/(1-f^2)^(3/2)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f     = root_2.first;
        const double f_x   = root_11.second;
        const double f_y   = root_12.second;
        const double f_xy  = root_2.second;
        const double sqrtf = 1./std::sqrt (std::abs (1.-f*f)+TOLERANCE);
        if(errno)
        {
            errors.add_error(errno, "acos");
        }

        const double sqrt3f = sqrtf/(std::abs (1.-f*f)+TOLERANCE);
        if(errno)
        {
            errors.add_error(errno, "acos");
        }

        const double df      = -sqrtf*f_xy-f*sqrt3f*f_x*f_y;
        return std::make_pair (std::asin (f), df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("acos (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class acos_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "acos";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new acos_FunctionNode (expr);
    }
};

// Arc Tangent
//--------------------------------------------------------------------------
class atan_FunctionNode : public FunctionNode
{
public:
    atan_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = atan (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "atan");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=atan (f) ==> dh/df=1/(1+f^2)df/dx
         */
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = atan (root.first);
        double dvalue = root.second/(1.+root.first*root.first);

        if(errno)
        {
            errors.add_error(errno, "atan");
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=atan (f) ==>d^2h/dxdy = d (1/(1+f^2)df/dy)/dx
         * = 1/(1+f^2)*d^2f/dxdy - 2f/(1+f^2)^2*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f     = root_2.first;
        const double f_x   = root_11.second;
        const double f_y   = root_12.second;
        const double f_xy  = root_2.second;
        const double f_1   = 1./(1.+f*f);
        if(errno)
        {
            errors.add_error(errno, "atan");
        }

        const double f_2 = f_1*f_1;
        if(errno)
        {
            errors.add_error(errno, "atan");
        }

        const double df      = f_1*f_xy-2.*f*f_2*f_x*f_y;
        return std::make_pair (std::atan (f), df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("atan (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class atan_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "atan";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new atan_FunctionNode (expr);
    }
};

// Arc Tangent 2: atan2(y, x)
//--------------------------------------------------------------------------
class atan2_FunctionNode : public FunctionNode
{
public:
    atan2_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double result = std::atan2(m_nodes[0]->Evaluate (errors), m_nodes[1]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "atan2");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'atan2' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'atan2' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("atan2 (" + m_nodes[0]->get_renamed_expression (var_map_name) + ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class atan2_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "atan2";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new atan2_FunctionNode (expr);
    }
};

// Log
// means log10
//--------------------------------------------------------------------------
class log_FunctionNode : public FunctionNode
{
public:
    log_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;
        double tmp = m_nodes[0]->Evaluate (errors);

        double result = std::log10(tmp);

        if(errno)
        {
            errors.add_error(errno, "log");
        }

        if(tmp<TOLERANCE)
        {
            result = -INF_VALUE;
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=log10(f) ==>dh/dx = 1/(log10*f)*df/dx
         */
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);

        double value = std::log10(root.first);
        if(errno)
        {
            errors.add_error(errno, "log");
        }
        double dvalue = 0.434294481903252*root.second/(std::abs (root.first)+TOLERANCE); // 1/ln (10) = 0.434294481903252
        if(errno)
        {
            errors.add_error(errno, "log");
        }

        if(root.first<TOLERANCE)
        {
            value = -INF_VALUE;
            dvalue = INF_VALUE*sign (root.second);
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=log(f) ==>d^2h/dxdy = d (1/log10*f^(-1)df/dy)/dx
         * = 1/log10*f^(-1)*d^2f/dxdy - 1/log10*f^(-2)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f     = root_2.first;
        const double f_x   = root_11.second;
        const double f_y   = root_12.second;
        const double f_xy  = root_2.second;
        const double f_1   = 1./(f+TOLERANCE*sign (f));
        const double f_2   =  f_1*f_1;

        double df    = 0.434294481903252*(f_1*f_xy-f_2*f_x*f_y);
        double logf  = std::log10(f);
        if(errno)
        {
            errors.add_error(errno, "log");
        }

        if(f<TOLERANCE)
        {
            logf = -INF_VALUE;
            df   = -INF_VALUE*sign (f_xy);
        }
        return std::make_pair (logf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("log (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class log_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "log";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new log_FunctionNode (expr);
    }
};

// Ln
//--------------------------------------------------------------------------
class ln_FunctionNode : public FunctionNode
{
public:
    ln_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;
        double tmp = m_nodes[0]->Evaluate (errors);

        double result = std::log(tmp);

        if(errno)
        {
            errors.add_error(errno, "ln");
        }

        if(tmp<TOLERANCE)
        {
            result = -INF_VALUE;
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double value = log(root.first);
        double dvalue = root.second/(root.first+1.e-20);

        if(errno)
        {
            errors.add_error(errno, "ln");
        }

        if(root.first<TOLERANCE)
        {
            value = -INF_VALUE;
            dvalue = INF_VALUE*sign (root.second);
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=ln (f) ==>d^2h/dxdy = d (f^(-1)df/dy)/dx
         * = f^(-1)*d^2f/dxdy - f^(-2)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f     = root_2.first;
        const double f_x   = root_11.second;
        const double f_y   = root_12.second;
        const double f_xy  = root_2.second;
        const double f_1   = 1./(f+TOLERANCE*sign (f));
        const double f_2   =  f_1*f_1;

        double df    = f_1*f_xy-f_2*f_x*f_y;
        double logf  = std::log(f);
        if(errno)
        {
            errors.add_error(errno, "log");
        }

        if(f<TOLERANCE)
        {
            logf = -INF_VALUE;
            df   = -INF_VALUE*sign (f_xy);
        }
        return std::make_pair (logf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("ln (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class ln_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "ln";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new ln_FunctionNode (expr);
    }
};

// Exp
//--------------------------------------------------------------------------
class exp_FunctionNode : public FunctionNode
{
public:
    exp_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        const double exponent = m_nodes[0]->Evaluate (errors);

        double result = 0.;

        if (exponent < std::numeric_limits<double>::min_exponent*log(2.))
        {
            result = 0.;
        }
        else if (exponent > std::numeric_limits<double>::max_exponent*log(2.))
        {
            result = INF_VALUE;
        }
        else
        {
            result = std::exp(exponent);

            if(errno)
            {
                errors.add_error(errno, "exp");
            }
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=exp(f) ==> dh/dx = exp(f)df/dx
         */
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        const double exponent = root.first;
        double value = 0.;

        if(exponent < std::numeric_limits<double>::min_exponent*log(2.))
        {
            value = 0.;
        }
        else if(exponent > std::numeric_limits<double>::max_exponent*log(2.))
        {
            value = INF_VALUE;
        }
        else
        {
            value = std::exp(exponent);

            if(errno)
            {
                errors.add_error(errno, "exp");
            }
        }

        double dvalue = root.second*value;

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=exp(f) ==>d^2h/dxdy = d (exp(f)df/dy)/dx
         * = exp(f)*(d^2f/dxdy + df/dx*df/dy)
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        double expf = 0.;//

        if(f<std::numeric_limits<double>::min_exponent*log(2.))
        {
            expf = 0.;
        }
        else if(f>std::numeric_limits<double>::max_exponent*log(2.))
        {
            expf = INF_VALUE;
        }
        else
        {
            expf = std::exp(f);
            if(errno)
            {
                errors.add_error(errno, "exp");
            }
        }

        const double df = expf*(f_xy+f_x*f_y);
        return std::make_pair (expf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("exp (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class exp_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "exp";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new exp_FunctionNode (expr);
    }
};

// Logn
//--------------------------------------------------------------------------
class logn_FunctionNode : public FunctionNode
{
public:
    logn_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        // Check for division by zero
        double tmp1 = m_nodes[1]->Evaluate (errors);
        double tmp = log(tmp1);

        if(errno)
        {
            errors.add_error(errno, "logn");
        }

        // Calculate result
        double tmp0=m_nodes[0]->Evaluate (errors);
        double result = log(tmp0) / (tmp+TOLERANCE*sign (tmp));

        if(errno)
        {
            errors.add_error(errno, "logn");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        // Check for division by zero
        double tmp1 = m_nodes[1]->Evaluate (errors);
        double tmp = log(tmp1);
        tmp+=TOLERANCE*sign (tmp);
        if(errno)
        {
            errors.add_error(errno, "logn");
        }

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);

        double value = log(root.first)/tmp;
        double dvalue = root.second/(std::abs (root.first)+TOLERANCE)/tmp;

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        /*
         * h=logn (f) ==>d^2h/dxdy = d (1/logn*f^(-1)df/dy)/dx
         * = 1/logn*f^(-1)*d^2f/dxdy - 1/logn*f^(-2)*df/dx*df/dy
         */
        const double exponent = m_nodes[1]->Evaluate (errors);
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f     = root_2.first;
        const double f_x   = root_11.second;
        const double f_y   = root_12.second;
        const double f_xy  = root_2.second;
        const double f_1   = 1./(f+TOLERANCE*sign (f));
        const double f_2   =  f_1*f_1;

        const double logn = log(exponent);
        if(errno)
        {
            errors.add_error(errno, "logn");
        }
        const double logn_1 = 1./(logn+TOLERANCE*sign (logn));

        double df    = logn_1*(f_1*f_xy-f_2*f_x*f_y);
        double logf  = std::log(f)*logn_1;
        if(errno)
        {
            errors.add_error(errno, "log");
        }

        if(f<TOLERANCE)
        {
            logf = -INF_VALUE;
            df   = -INF_VALUE*sign (f_xy);
        }
        return std::make_pair (logf, df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("logn (" + m_nodes[0]->get_renamed_expression (var_map_name) + ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class logn_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "logn";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new logn_FunctionNode (expr);
    }
};

// Ceil
//--------------------------------------------------------------------------
class ceil_FunctionNode : public FunctionNode
{
public:
    ceil_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        return ceil(m_nodes[0]->Evaluate (errors));
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'ceil' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'ceil' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("ceil (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class ceil_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "ceil";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new ceil_FunctionNode (expr);
    }
};

// Floor
//--------------------------------------------------------------------------
class floor_FunctionNode : public FunctionNode
{
public:
    floor_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        return floor(m_nodes[0]->Evaluate (errors));
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'floor' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'floor' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("floor (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class floor_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "floor";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new floor_FunctionNode (expr);
    }
};

// Rand
//--------------------------------------------------------------------------

// Get next random (0,1)
inline double NextRandom(double *seed)
{
    long a = (long)(*seed) * 214013L + 2531011L;
    *seed = (double)a;

    return (double)((a >> 16) & 0x7FFF) / 32767.0;
}

class rand_FunctionNode : public FunctionNode
{
public:
    rand_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (0, 0, 1, 1);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        (void)errors;
        return NextRandom(m_refs[0]);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rand' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rand' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("rand (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class rand_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "rand";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new rand_FunctionNode (expr);
    }
};

// Random
//--------------------------------------------------------------------------
class random_FunctionNode : public FunctionNode
{
public:
    random_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 1, 1);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double a = m_nodes[0]->Evaluate (errors);
        double b = m_nodes[1]->Evaluate (errors);

        return NextRandom(m_refs[0]) * (b - a) + a;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'random' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'random' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("random (" + m_nodes[0]->get_renamed_expression (var_map_name) + ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class random_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "random";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new random_FunctionNode (expr);
    }
};

// Randomize
//--------------------------------------------------------------------------
class randomize_FunctionNode : public FunctionNode
{
public:
    randomize_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (0, 0, 1, 1);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        (void)errors;
        static long curcall = 1;

        *m_refs[0] = (double)(((clock() + 1024u + curcall) * time (NULL)) % 2176971487u);
        curcall++;

        return 0.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'randomize' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'randomize' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("randomize (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class randomize_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "randomize";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new randomize_FunctionNode (expr);
    }
};

// Radians to degrees
//--------------------------------------------------------------------------
class deg_FunctionNode : public FunctionNode
{
public:
    deg_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        return (m_nodes[0]->Evaluate (errors) * 180.0) / M_PI;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double tmp = 180./M_PI;
        double value = root.first*tmp;
        double dvalue = root.second*tmp;

        if(errno)
            SHOW_FATAL_ERROR("deg");

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'deg' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("deg (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class deg_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "deg";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new deg_FunctionNode (expr);
    }
};

// Degrees to radians
//--------------------------------------------------------------------------
class rad_FunctionNode : public FunctionNode
{
public:
    rad_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        return (m_nodes[0]->Evaluate (errors) * M_PI) / 180.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;

        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);
        double tmp = M_PI/180.;
        double value = root.first*tmp;
        double dvalue = root.second*tmp;

        if(errno)
            SHOW_FATAL_ERROR("rad");

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rad' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("rad (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class rad_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "rad";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new rad_FunctionNode (expr);
    }
};

// Rectangular to polar: rect2pol(x, y, &distance, &angle)
//--------------------------------------------------------------------------
class rect2pol_FunctionNode : public FunctionNode
{
public:
    rect2pol_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 2, 2);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double x = m_nodes[0]->Evaluate (errors);
        double y = m_nodes[1]->Evaluate (errors);

        double d = sqrt (x * x + y * y);
        double a = atan2(y, x);

        if(errno)
            SHOW_FATAL_ERROR("rect2pol");

        *m_refs[0] = d;
        if(a < 0.0)
            *m_refs[1] = a + (2.0 * M_PI);
        else
            *m_refs[1] = a;

        return d;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rect2pol' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rect2pol' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("rect2pol (" + m_nodes[0]->get_renamed_expression (var_map_name) + ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class rect2pol_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "rect2pol";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new rect2pol_FunctionNode (expr);
    }
};

// Polar to rectangular: pol2rect (distance, angle, &x, &y)
//--------------------------------------------------------------------------
class pol2rect_FunctionNode : public FunctionNode
{
public:
    pol2rect_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 2, 2);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        double d = m_nodes[0]->Evaluate (errors);
        double a = m_nodes[1]->Evaluate (errors);

        double x = d * cos (a);
        double y = d * sin (a);

        if(errno)
            SHOW_FATAL_ERROR("pol2rect");

        *m_refs[0] = x;
        *m_refs[1] = y;

        return x;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'pol2rect' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'pol2rect' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("pol2rect (" + m_nodes[0]->get_renamed_expression (var_map_name) + ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class pol2rect_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "pol2rect";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new pol2rect_FunctionNode (expr);
    }
};

// If
//--------------------------------------------------------------------------
class if_FunctionNode : public FunctionNode
{
public:
    if_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (3, 3, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double c = m_nodes[0]->Evaluate (errors);

        if(equal_eps (c, 0.0))
            return m_nodes[2]->Evaluate (errors);
        else
            return m_nodes[1]->Evaluate (errors);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'if' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'if' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("if (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[2]->get_renamed_expression (var_map_name) + ")");
    }
};

class if_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "if";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new if_FunctionNode (expr);
    }
};

// Select
//--------------------------------------------------------------------------
class select_FunctionNode : public FunctionNode
{
public:
    select_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (3, 4, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double c = m_nodes[0]->Evaluate (errors);

        if(equal_eps (c, 0.0))
            return m_nodes[2]->Evaluate (errors);
        else if(c < 0.0)
            return m_nodes[1]->Evaluate (errors);
        else
        {
            if(m_nodes.size() == 3)
                return m_nodes[2]->Evaluate (errors);
            else
                return m_nodes[3]->Evaluate (errors);
        }
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'select' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'select' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("select (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class select_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "select";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new select_FunctionNode (expr);
    }
};

// Equal
//--------------------------------------------------------------------------
class equal_FunctionNode : public FunctionNode
{
public:
    equal_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        if(equal_eps (m_nodes[0]->Evaluate (errors), m_nodes[1]->Evaluate (errors)))
            return 1.0;
        else
            return 0.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'equal' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'equal' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("equal (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class equal_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "equal";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new equal_FunctionNode (expr);
    }
};

// Above
//--------------------------------------------------------------------------
class above_FunctionNode : public FunctionNode
{
public:
    above_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        if(m_nodes[0]->Evaluate (errors) > m_nodes[1]->Evaluate (errors))
            return 1.0;
        else
            return 0.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'above' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'above' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("above (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class above_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "above";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new above_FunctionNode (expr);
    }
};

// Below
//--------------------------------------------------------------------------
class below_FunctionNode : public FunctionNode
{
public:
    below_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        if(m_nodes[0]->Evaluate (errors) < m_nodes[1]->Evaluate (errors))
            return 1.0;
        else
            return 0.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'below' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'below' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("below (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class below_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "below";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new below_FunctionNode (expr);
    }
};

// Clip
//--------------------------------------------------------------------------
class clip_FunctionNode : public FunctionNode
{
public:
    clip_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (3, 3, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double v = m_nodes[0]->Evaluate (errors);
        double a = m_nodes[1]->Evaluate (errors);
        double b = m_nodes[2]->Evaluate (errors);

        if(v < a)
            return a;
        else if(v > b)
            return b;
        else
            return v;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'clip' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'clip' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("clip (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[2]->get_renamed_expression (var_map_name) + ")");
    }
};

class clip_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "clip";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new clip_FunctionNode (expr);
    }
};

// Clamp
//--------------------------------------------------------------------------
class clamp_FunctionNode : public FunctionNode
{
public:
    clamp_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (3, 3, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double v = m_nodes[0]->Evaluate (errors);
        double a = m_nodes[1]->Evaluate (errors);
        double b = m_nodes[2]->Evaluate (errors);

        if(equal_eps (a, b))
            return a;
        else
        {
            double tmp = fmod (v - a, b - a);

            if(tmp < 0)
                return tmp + b;
            else
                return tmp + a;
        }
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'clamp' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'clamp' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("clamp (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[2]->get_renamed_expression (var_map_name) + ")");
    }
};

class clamp_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "clamp";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new clamp_FunctionNode (expr);
    }
};

// Rescale
//--------------------------------------------------------------------------
class rescale_FunctionNode : public FunctionNode
{
public:
    rescale_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (5, 5, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double pnt = m_nodes[0]->Evaluate (errors);
        double o1 = m_nodes[1]->Evaluate (errors);
        double o2 = m_nodes[2]->Evaluate (errors);
        double n1 = m_nodes[3]->Evaluate (errors);
        double n2 = m_nodes[4]->Evaluate (errors);

        double odiff = o2 - o1;
        if(equal_eps (odiff, 0.0))
            return n1;
        else
        {
            return (pnt - o1) * (n2 - n1) / odiff + n1;
        }
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rescale' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'rescale' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("rescale (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class rescale_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "rescale";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new rescale_FunctionNode (expr);
    }
};

// Poly: poly (x, c3, c2, c1, c0): c3*x^3 + c2*x^2 + c1*x + c0
//--------------------------------------------------------------------------
class poly_FunctionNode : public FunctionNode
{
public:
    poly_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, -1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        double total = 0.0;
        double curpow;

        std::vector<Node*>::size_type pos, count;
        count = m_nodes.size();

        curpow = (double)(count - 2);

        // Value of x
        double x = m_nodes[0]->Evaluate (errors);

        errno = 0;

        for (pos = 1; pos < count; pos++)
        {
            total += (m_nodes[pos]->Evaluate (errors) * pow(x, curpow));
            curpow -= 1.0;
        }

        if(errno)
            SHOW_FATAL_ERROR("poly");

        return total;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        double value = 0.0;
        double dvalue = 0.;
        double curpow;

        std::vector<Node*>::size_type pos, count;
        count = m_nodes.size();

        curpow = (double)(count - 2);

        // Value of x
        std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);

        errno = 0;

        for (pos = 1; pos < count; pos++)
        {
            double coe = m_nodes[pos]->Evaluate (errors);
            value += ( coe * pow(root.first, curpow));
            if(pos!=count-1)
            {
                dvalue += coe * curpow * pow(root.first, curpow - 1.)*root.second;
            }
            curpow -= 1.0;
        }

        if(errno)
            SHOW_FATAL_ERROR("poly");

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'poly' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("poly (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class poly_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "poly";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new poly_FunctionNode (expr);
    }
};

// And
//--------------------------------------------------------------------------
class and_FunctionNode : public FunctionNode
{
public:
    and_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        if(equal_eps (m_nodes[0]->Evaluate (errors), 0.0) || equal_eps (m_nodes[1]->Evaluate (errors), 0.0))
            return 0.0;
        else
            return 1.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'and' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'and' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("and (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class and_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "and";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new and_FunctionNode (expr);
    }
};

// Or
//--------------------------------------------------------------------------
class or_FunctionNode : public FunctionNode
{
public:
    or_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (2, 2, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        if(equal_eps (m_nodes[0]->Evaluate (errors), 0.0) && equal_eps (m_nodes[1]->Evaluate (errors), 0.0))
            return 0.0;
        else
            return 1.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'or' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'or' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("or (" + m_nodes[0]->get_renamed_expression (var_map_name) +
                ", " + m_nodes[1]->get_renamed_expression (var_map_name) + ")");
    }
};

class or_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "or";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new or_FunctionNode (expr);
    }
};

// Not
//--------------------------------------------------------------------------
class not_FunctionNode : public FunctionNode
{
public:
    not_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        if(equal_eps (m_nodes[0]->Evaluate (errors), 0.0))
            return 1.0;
        else
            return 0.0;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'not' function is not supported when derivativing it");
        (void)var_symbol;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        SHOW_FATAL_ERROR("'not' function is not supported when derivativing it");
        (void)var_symbol_1;
        (void)var_symbol_2;
        (void)errors;
        return std::make_pair (0., 0.);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("not (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class not_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "not";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new not_FunctionNode (expr);
    }
};


// erf
//--------------------------------------------------------------------------
class erf_FunctionNode : public FunctionNode
{
public:
    erf_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        const double result = erf (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "erf");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;
        const std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);

        const double value = erf (root.first);
        if(errno)
        {
            errors.add_error(errno, "erf");
        }

        /*
         * h=erf (f) ==> dh/dx = 2/sqrt(pi)*exp(-f^2)df/dx
         */
        double dvalue = root.second * 1.1284 * exp (-(root.first * root.first));
        if(errno)
        {
            errors.add_error(errno, "erf");
        }
        if(errno == ERANGE)
        {
            dvalue = INF_VALUE;
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=erf (f) ==> d^2h/dxdy = 2/sqrt(pi)*exp(-f^2)*d^2f/dxdy - 4/sqrt(pi)*f*exp(-f^2)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double df   = 1.1284 * exp(-f*f)*f_xy - 2.2568*f*exp(-f*f)*f_x*f_y;
        return std::make_pair (erf (f), df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("erf (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class erf_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "erf";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new erf_FunctionNode (expr);
    }
};


// erfc
//--------------------------------------------------------------------------
class erfc_FunctionNode : public FunctionNode
{
public:
    erfc_FunctionNode (Expression *expr) : FunctionNode (expr)
    {
        SetArgumentCount (1, 1, 0, 0);
    }

    double DoEvaluate (ErrorInfo& errors) const
    {
        errno = 0;

        const double result = erfc (m_nodes[0]->Evaluate (errors));

        if(errno)
        {
            errors.add_error(errno, "erf");
        }

        return result;
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol, ErrorInfo& errors) const
    {
        errno = 0;
        const std::pair<double, double> root = m_nodes[0]->Evaluate (var_symbol, errors);

        const double value = erfc (root.first);
        if(errno)
        {
            errors.add_error(errno, "erf");
        }

        /*
         * h=erf (f) ==> dh/dx = 2/sqrt(pi)*exp(-f^2)df/dx
         */
        double dvalue = -(root.second * 1.1284 * exp (-(root.first * root.first)));
        if(errno)
        {
            errors.add_error(errno, "erf");
        }
        if(errno == ERANGE)
        {
            dvalue = INF_VALUE;
        }

        return std::make_pair (value, dvalue);
    }

    std::pair<double, double> DoEvaluate (const std::string& var_symbol_1, const std::string& var_symbol_2, ErrorInfo& errors) const
    {
        errno = 0;

        /*
         * h=erf (f) ==> d^2h/dxdy = 2/sqrt(pi)*exp(-f^2)*d^2f/dxdy - 4/sqrt(pi)*f*exp(-f^2)*df/dx*df/dy
         */
        const std::pair<double, double> root_2  = m_nodes[0]->Evaluate (var_symbol_1, var_symbol_2, errors);
        const std::pair<double, double> root_11 = m_nodes[0]->Evaluate (var_symbol_1, errors);
        const std::pair<double, double> root_12 = m_nodes[0]->Evaluate (var_symbol_2, errors);
        const double f    = root_2.first;
        const double f_x  = root_11.second;
        const double f_y  = root_12.second;
        const double f_xy = root_2.second;
        const double df   = -1.1284 * exp(-f*f)*f_xy + 2.2568*f*exp(-f*f)*f_x*f_y;
        return std::make_pair (erfc(f), df);
    }

    std::string
    get_renamed_expression (
        const std::string& var_map_name) const
    {
        return ("erfc (" + m_nodes[0]->get_renamed_expression (var_map_name) + ")");
    }
};

class erfc_FunctionFactory : public FunctionFactory
{
public:
    std::string GetName() const
    {
        return "erfc";
    }

    FunctionNode *DoCreate (Expression *expr)
    {
        return new erfc_FunctionNode (expr);
    }
};

// Initialize default functions
void FunctionList::AddDefaultFunctions()
{
#define ADDFUNCTION(name) \
        std::auto_ptr<FunctionFactory> name ## _func(new name ## _FunctionFactory()); \
        m_functions.push_back(name ## _func.get()); \
        name ## _func.release();

    ADDFUNCTION (abs);
    ADDFUNCTION (mod);

    ADDFUNCTION (ipart);
    ADDFUNCTION (fpart);

    ADDFUNCTION (min);
    ADDFUNCTION (max);
    ADDFUNCTION (sqrt);

    ADDFUNCTION (sin);
    ADDFUNCTION (cos);
    ADDFUNCTION (tan);

    ADDFUNCTION (sinh);
    ADDFUNCTION (cosh);
    ADDFUNCTION (tanh);

    ADDFUNCTION (asin);
    ADDFUNCTION (acos);
    ADDFUNCTION (atan);
    ADDFUNCTION (atan2);

    ADDFUNCTION (log);
    ADDFUNCTION (ln);
    ADDFUNCTION (exp);
    ADDFUNCTION (logn);

    ADDFUNCTION (ceil);
    ADDFUNCTION (floor);

    ADDFUNCTION (rand);
    ADDFUNCTION (random);
    ADDFUNCTION (randomize);

    ADDFUNCTION (deg);
    ADDFUNCTION (rad);

    ADDFUNCTION (rect2pol);
    ADDFUNCTION (pol2rect);

    ADDFUNCTION (if); // Preprocess will take care of this beforehand
    ADDFUNCTION (select);

    ADDFUNCTION (equal);
    ADDFUNCTION (above);
    ADDFUNCTION (below);

    ADDFUNCTION (clip);
    ADDFUNCTION (clamp);
    ADDFUNCTION (rescale);

    ADDFUNCTION (poly);

    ADDFUNCTION (and);
    ADDFUNCTION (or);
    ADDFUNCTION (not);

    ADDFUNCTION (erf);
    ADDFUNCTION (erfc);
}

}
